function [ code_control,places_grid,gridmap,edges,places_gecon ] = brazil_grid( code_control,country_bounds,places_sedac,...
                                                                                   relief_etopo,places_gecon)                                                                            
% this file REPLACES create_grid_map
datafolder_LAC = ['../../data/LAC/'];
datafolder_LAC_country = ['../../data/LAC/BR/'];

% load in a shapfile of the "meso regions" of Brazil
filepath = [datafolder_LAC_country, 'br_mesorregioes/BR_Mesorregioes_2019.shp'];
brazil_regions = shaperead(filepath); 

% simplify each polygon 
for row=1:length(brazil_regions)
    [ row length(brazil_regions) ]
    original = length( brazil_regions(row).Y );
    
    brY = brazil_regions(row).Y;
    brX = brazil_regions(row).X;
    [brY1, brX1, cerr, tol] = reducem(brY', brX');
    brazil_regions(row).Y = brY1';
    brazil_regions(row).X = brX1';
    
    new = length(  brazil_regions(row).Y );
    [ original new ];
end

% overlay the G-ECON data to compute value added by mesoregion
load([datafolder_LAC_country, 'popplaces_sedac.mat'])
load([datafolder_LAC_country, 'relief_etopo.mat'])
load([datafolder_LAC_country, 'country_bounds.mat'])

%% population, SEDAC 
places_X = cell2mat({places_sedac.X}');
places_Y = cell2mat({places_sedac.Y}');
population = cell2mat({places_sedac.population}');

for row=1:length(brazil_regions)
    [ row length(brazil_regions) ]
    in_region = inpolygon(places_X, places_Y, brazil_regions(row).X, brazil_regions(row).Y);
    brazil_regions(row).newPop = sum(in_region.*population);
end

%% import official statistics from stats webpage
official_stats = table2struct(readtable([datafolder_LAC, 'BR/official_stats_meso.csv']));

% shapefile has two extra, somehow.. 4377 and 4388
k = 1;
j = length(brazil_regions);
while k <= j
   if strcmp(brazil_regions(k).CD_MESO,'4377') || strcmp(brazil_regions(k).CD_MESO,'4388')
        brazil_regions(k) = [];
   else 
        k = k+1;
   end
   j = length(brazil_regions);
end

%% COMBINING THESE INTO THE NEW GRID
gridmap = [];
for i=1:length(brazil_regions)
        gridmap(i).Geometry = 'Polygon';
        gridmap(i).X = brazil_regions(i).X;
        gridmap(i).Y = brazil_regions(i).Y;
        gridmap(i).counter = i;
        gridmap(i).mesoregion = brazil_regions(i).CD_MESO;
        
        gridmap(i).population = brazil_regions(i).newPop;
        gridmap(i).pop2 =   official_stats(i).pop;
        gridmap(i).income = official_stats(i).gdp;
        
        if gridmap(i).population == 0
            gridmap(i).income = 0;
        end
        
        gridmap(i).y = gridmap(i).income/gridmap(i).population;  %income per capita
        if isnan(gridmap(i).y)
            gridmap(i).y = 0;
            gridmap(i).y = cast(gridmap(i).y,'single');
        end
        
        i=i+1;
end

%{ 
    forest extent -- copying format from etopo
if code_control.brazil_forest == 1
    file  = [ datafolder_LAC_country,'gfc.mat'];
    load( file,'gfc' );    

    tree_cover = cell2mat( {gfc.tree_cover} );
    places_X_gfc = cell2mat({gfc.centerX}');
    places_Y_gfc = cell2mat({gfc.centerY}');

    % these two regions are islands once we eliminate unbuildable areas 
    exclude = [2, 9];
    for j=1:length( gridmap )

            % roughness
            in_j = inpolygon( places_X_gfc, places_Y_gfc, gridmap(j).X,gridmap(j).Y  ); % indicate whether each 0.25 degree cell in Etopo belongs in the cell

            gridmap(j).tree_cover = mean( tree_cover( in_j ) ); 
            gridmap(j).build = (gridmap(j).tree_cover < 85 || isnan( gridmap(j).tree_cover )) ...
                && ~ismember(j, exclude);

    end
else
%}
    for j=1:length( gridmap )
        gridmap(j).build = 1;
    end
% end

% eliminate mesoregions that don't meet the criteria
k = length(gridmap);
j = 1;
while j <= k
   if gridmap(j).build == 0
       gridmap(j) = [];
   else
        j = j + 1;
   end
   k = length(gridmap);
end

%% ruggedness, topographic variables
places_X_topo = cell2mat({relief_etopo.centerX}');
places_Y_topo = cell2mat({relief_etopo.centerY}');

av_alt = cell2mat( {relief_etopo.av_alt} );
sd_alt = cell2mat( {relief_etopo.sd_alt} );
rugged = cell2mat( {relief_etopo.rugged} );

for j=1:length( gridmap )

        % roughness
        in_j = inpolygon( places_X_topo, places_Y_topo, gridmap(j).X,gridmap(j).Y  ); % indicate whether each 0.25 degree cell in Etopo belongs in the cell
        
        gridmap(j).av_alt = mean( av_alt( in_j ) );  
        gridmap(j).sd_alt = mean( sd_alt( in_j ) );  
        gridmap(j).rugged = mean( rugged( in_j ) ); 

end

% verify population from G-econ matches with population from SEDAC
plot( cell2mat( {gridmap.population} ) ,cell2mat( {gridmap.pop2} ),'or' )
corr( cell2mat( {gridmap.population} )',cell2mat( {gridmap.pop2} )' )

% determine centroids
% X and Y coordinates of each city
places_X = cell2mat({places_sedac.X}');
places_Y = cell2mat({places_sedac.Y}');

% population vector
pop_mat = cell2mat( {places_sedac.population}' ); % vector with city population

% recompute areas and population
keep = ones( length(gridmap), 1);
for j=1:length( gridmap )
        [ j ]
    if keep(j)==1
        
        % total area
        in_j = inpolygon( places_X, places_Y, gridmap(j).X,gridmap(j).Y  ); % indicator that each city is in cell j
        gridmap(j).area = sum( areaint( gridmap(j).Y,gridmap(j).X ) );
        
        if code_control.pick_largest   % pick largest city in each cell

            [ ~,j_largest ] = max( in_j.*pop_mat );  % index of largest city in cell j
            gridmap( j ).coordinates = [ places_X( j_largest ) places_Y( j_largest ) ];

        else  % pick centroids

            gridmap( j ).coordinates = [ ( inpolygon( places_X, places_Y, gridmap(j).X,gridmap(j).Y  ).*places_X )'*pop_mat/gridmap(j).population ...
                                         ( inpolygon( places_X, places_Y, gridmap(j).X,gridmap(j).Y  ).*places_Y )'*pop_mat/gridmap(j).population  ];           

            % centroid can be outside the country for non-convex regions, if so pick largest city
            if ~inpolygon( gridmap(j).coordinates(1),gridmap(j).coordinates(2),country_bounds.X,country_bounds.Y )          
                    [ ~,j_largest ] = max( in_j.*pop_mat );  % index of largest city in cell j
                    gridmap( j ).coordinates = [ places_X( j_largest ) places_Y( j_largest ) ];   
            end
        end
    end
    
end

% define a neighbor as a node that is within X miles 
centroidX = zeros(length(gridmap), 1);
centroidY = zeros(length(gridmap), 1);
for j=1:length(gridmap)
    centroidX(j) = gridmap( j ).coordinates(1);
    centroidY(j) = gridmap( j ).coordinates(2);
end

% determine neighbors 
% first, compute a distance matrix
dist_mat = ones(length(gridmap), length(gridmap));
for j = 1:length(gridmap)
    for k=(j+1):length(gridmap)
        [ j k ]
        dist_mat(j, k) = deg2km(distance(centroidX(j),centroidY(j),centroidX(k),centroidY(k)));
        dist_mat(k, j) = dist_mat(j, k);
    end
end

% or neighboring cells (shared border)
border = zeros(length(gridmap));
struct2poly = @(s) polyshape(s.X,s.Y); % define my struct -> polyshape operation
p = arrayfun(struct2poly,gridmap);
pfat = polybuffer(p,1e-4);

for k = 1:length(gridmap)
    for j = (k+1):length(gridmap) % Only need j for ks we haven't checked yet
        [ k j ]
        border(j,k) = area(intersect(pfat(j),pfat(k))) > 1e-6;
    end
end

for k = 1:length(gridmap)
    gridmap(k).neighbors = [];
    for j = (k+1):length(gridmap) 
        [ j k ]
        if border(j,k) == 1 && (isempty(gridmap(k).neighbors)==1 ...
                || max(ismember(gridmap(k).neighbors, j))==0 )
            if dist_mat(k, j) < 800
                gridmap(k).neighbors = [gridmap(k).neighbors j];
            end            
        end
    end
end

% make a dataset of the edges just for the sake of understanding
% and seeing which looks better
% how many edges on each graph?
% vectors with the origin and destination nodes
gridmap_s = [];
gridmap_t = [];
for j=1:length(gridmap)
    gridmap(j).counter = j; % reset the counter bc of the deleted rows
    gridmap_s = [ gridmap_s;gridmap(j).counter*ones( length( gridmap(j).neighbors ),1 ) ];
    gridmap_t = [ gridmap_t;gridmap(j).neighbors' ];
end
temp = [ gridmap_s gridmap_t];

% keep unique horizontal-vertical edges
for i=1:length( temp )    
    temp(i,:) = [ min( temp(i,:) ) max( temp(i,:) ) ];    
end    
[ unique_edges,~,~ ] = unique( temp,'rows');  % we will use unique_edges to build the vertical-horizontal graph

% population
temp =  cell2mat({ gridmap.population })';
quantiles_pop_grid = quantile( temp( temp>0 ),code_control.N_cat );  % assign quantiles conditional on positive population in the square
for i=1:length(gridmap)
    if temp(i)>0
        [ gridmap(i).quantiles ] = sum( gridmap(i).population>quantiles_pop_grid )+1;
    elseif temp(i)==0
        [ gridmap(i).quantiles ] = 1; % changed from 0
    end
end

% income
temp =  cell2mat({ gridmap.y })';
quantiles_income_grid = quantile( temp( temp>0 ),code_control.N_cat );  % assign quantiles conditional on positive population in the square
for i=1:length(gridmap)
    if temp(i)>0
        [ gridmap(i).quantiles_inc ] = sum( gridmap(i).y>quantiles_income_grid )+1;
    elseif temp(i)==0
        [ gridmap(i).quantiles_inc ] = 1; % changed from 0
    end
end   

% assign altitude and ruggedness
temp1 =  cell2mat({ gridmap.av_alt })';
temp2 =  cell2mat({ gridmap.rugged })';
for i=1:length(gridmap)
    [ gridmap(i).quantiles_av_alt ] = sum( gridmap(i).av_alt>quantile( temp1,code_control.N_cat ) )+1;
    [ gridmap(i).quantiles_rugged ] = sum( gridmap(i).rugged>quantile( temp2,code_control.N_cat ) )+1;
end

% winsorize
ytemp = cell2mat( {gridmap.y} );
if length(gridmap)>20
    maxquantile = max( quantile(ytemp,20) );
elseif length(gridmap)>10
    maxquantile = max( quantile(ytemp,10) );
elseif length(gridmap)<=10
    maxquantile = max( quantile(ytemp,5) );
end

if code_control.winsorize == 1
    for i=1:length(gridmap)

        if gridmap(i).y>maxquantile
            gridmap(i).y = 0.99*maxquantile;
            gridmap(i).income = gridmap(i).y*gridmap(i).population;
        end 

    end
end

% create places
places_grid = [];
for j=1:length(gridmap)
    
    places_grid(j).Geometry = 'Point';
    places_grid(j).X = gridmap(j).coordinates(1);
    places_grid(j).Y = gridmap(j).coordinates(2);
    places_grid(j).population = gridmap(j).population;
    places_grid(j).income = gridmap(j).income;
    places_grid(j).y = gridmap(j).y;
    places_grid(j).quantiles = gridmap(j).quantiles;  %pop_quantile
    places_grid(j).neighbors = gridmap(j).neighbors;
    places_grid(j).av_alt = gridmap(j).av_alt;
    places_grid(j).sd_alt = gridmap(j).sd_alt;
    places_grid(j).rugged = gridmap(j).rugged;
    places_grid(j).build = gridmap(j).build;
    
end

edges.unique_edges = unique_edges;

%{ 
plotting 
addpath('../Toolbox/export_fig/');
graph_width=800;
graph_height=600;

% underlying grid 
close(figure(1))
fig = figure(1);
set(fig,'Position',[0 0 graph_width graph_height]);

mapshow(country_bounds, 'FaceColor','Black')
hold on;
mapshow( gridmap,...
           'FaceColor','Black',...
           'EdgeColor','White','LineWidth', 0.5 )
hold on;
mapshow( discretized_roads, 'Color',[ 1 0 0 ], ...
'LineWidth', 1  )
hold on;
N_cat = code_control.N_cat;
bright_nodes = 0.4; % higher=brighter
for i=N_cat+1:-1:1

    relsize = bright_nodes+(1-bright_nodes)*i/( N_cat+1 );
    
    mapshow( places_grid( cell2mat( {gridmap.quantiles} )==i ),...
        'Marker','o',...
        'MarkerFaceColor',relsize*[ 0 1 1 ],...
        'MarkerEdgeColor',[ 0 0 0 ],...
        'MarkerSize',6 )
end

margin = 0.5;
Xlim = [ min( cell2mat({country_bounds.X}) )-margin;max( cell2mat({country_bounds.X}) )+margin ];
Ylim = [ min( cell2mat({country_bounds.Y}) )-margin;max( cell2mat({country_bounds.Y}) )+margin ];
set( gca,'XTick',[],'YTick',[],'XLim',Xlim,'YLim',Ylim,'Box','on' )
export_fig([ path_papers_graph,country,'_grid.eps' ], '-depsc');

% population 
close(figure(2))
fig=figure(2);
set(fig,'Position',[0 0 graph_width graph_height]);

mapshow(country_bounds, 'FaceColor','Black')

% population
for i=N_cat+1:-1:1
    relsize = bright_nodes+(1-bright_nodes)*i/( N_cat+1 );
    mapshow( gridmap( cell2mat( {gridmap.quantiles} )==i ),...
        'FaceColor',relsize*[ 0 1 1 ] )
end

% centroids
for i=N_cat+1:-1:1

    relsize = bright_nodes+(1-bright_nodes)*i/( N_cat+1 );

    mapshow( places_grid( cell2mat( {gridmap.quantiles} )==i ),...
        'Marker','o',...
        'MarkerFaceColor',relsize*[ 0 1 1 ],...
        'MarkerEdgeColor',[ 0 0 0 ],...
        'MarkerSize',6 )
end

margin = 0.5;
Xlim = [ min( cell2mat({country_bounds.X}) )-margin;max( cell2mat({country_bounds.X}) )+margin ];
Ylim = [ min( cell2mat({country_bounds.Y}) )-margin;max( cell2mat({country_bounds.Y}) )+margin ];
set( gca,'XTick',[],'YTick',[],'XLim',Xlim,'YLim',Ylim,'Box','on' )
export_fig([ path_papers_graph,country,'_pop.eps' ], '-depsc');

% population 
close(figure(3))
fig=figure(3);
set(fig,'Position',[0 0 graph_width graph_height]);

mapshow(country_bounds, 'FaceColor','Black')

% population
for i=N_cat+1:-1:1
    relsize = bright_nodes+(1-bright_nodes)*i/( N_cat+1 );
    mapshow( gridmap( cell2mat( {gridmap.quantiles_inc} )==i ),...
        'FaceColor',relsize*[ 0 1 1 ] )
end

% centroids
for i=N_cat+1:-1:1

    relsize = bright_nodes+(1-bright_nodes)*i/( N_cat+1 );

    mapshow( places_grid( cell2mat( {gridmap.quantiles_inc} )==i ),...
        'Marker','o',...
        'MarkerFaceColor',relsize*[ 0 1 1 ],...
        'MarkerEdgeColor',[ 0 0 0 ],...
        'MarkerSize',6 )
end

margin = 0.5;
Xlim = [ min( cell2mat({country_bounds.X}) )-margin;max( cell2mat({country_bounds.X}) )+margin ];
Ylim = [ min( cell2mat({country_bounds.Y}) )-margin;max( cell2mat({country_bounds.Y}) )+margin ];
set( gca,'XTick',[],'YTick',[],'XLim',Xlim,'YLim',Ylim,'Box','on' )
export_fig([ path_papers_graph,country,'_income.eps' ], '-depsc');

%}






